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[57] ABSTRACT 

This is a method and associated apparatus for accu- 
rately and quickly estimating the amplitude, frequency 
and phase of a signal of interest. The method comprises 
the steps of, inputting the signal of interest; generating a 
reference signal with adjustable amplitude, frequency 
and phase at an output thereof; mixing the signal of 
interest with the reference signal and a signal 90° out of 
phase with the reference signal to provide a pair of 
quadrature sample signals comprising respectively a 
difference between the signal of interest and the refer- 
ence signal and a difference between the signal of inter- 
est and the signal 90° out of phase with the reference 
signal; using the pair of quadrature sample signals to 
compute estimates of the amplitude, frequency, and 
phase of an error signal comprising the difference be- 
tween the signal of interest and the reference signal 
employing a least squares estimation; adjusting the am- 
plitude, frequency, and phase of the reference signal 
from the numerically controlled oscillator in a manner 
which drives the error signal towards zero; and, output- 
ting the estimates of the amplitude, frequency, and 
phase of the error signal in combination with the refer- 
ence signal to produce a best estimate of the amplitude, 
frequency, and phase of the signal of interest. The pre- 
ferred method includes the step of providing the error 
signal as a real time confidence measure as to accuracy 
of the estimates wherein the closer the error signal is to 
zero, the higher the probability that the estimates are 
accurate. A matrix in the estimation algorithm provides 
an estimate of the variance of the estimation error. 

22 Claims, 2 Drawing Sheets 
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U.S. Pat. No. 4,533,918 of Vimot describes a device 


MODIFIED FAST FREQUENCY ACQUISITION working with radar equipment to obtain the location of 

VIA ADAPTIVE LEAST SQUARES ALGORITHM a target. This is achieved by a Kalman filter type of 

circuit that operates on the samples of range measure- 


ORIGIN ON THE INVENTION 5 

The invention described herein was made in the per- 
formance of work under a NASA contract, and is sub- 
ject to the provisions of Public Law 96-517 (35 USC 
202) in which the Contractor has elected not to retain 
title. 

This application is a continuation of application Ser. 

No. 523,692, filed May 15, 1990, now abandoned. 

TECHNICAL FIELD 

The invention relates to signal analysis apparatus and 15 
methods and, more particularly, a method of accurately 
and quickly estimating the amplitude, frequency and 
phase of a signal of interest. 

BACKGROUND ART 20 

The problem of estimating the parameters of a sinu- 
soidal signal has received considerable attention in the 
literature and in various prior art patents. Such a prob- 
lem arises in diverse engineering situations such as car- 25 
rier tracking for communications systems and the mea- 
surement of Doppler in position location, navigation, 
and radar systems. 

A variety of techniques have been proposed in the 
literature to solve such problems including, to mention 30 
a few, the application of the Fast Fourier Transform 
(FFT), one and two-dimensional Kalman filters based 
on a linearized model, a modified Kalman filter that 
results in a phase locked loop, and a digital phase locked 
loop derived on the basis of linear stochastic optimiza- 35 
tion. 

U.S. Pat. No. 4,686,532 pf McAulay describes a sys- 
tem for determining the strength of distributed scatter- 
ed (say P in number) that are possibly close to the sonar 
or radar relative to the array dimension of size I. As- 40 
suming that P scatterers rediate signals at some of the L 
Doppler frequencies known in advance, the method 
uses a least squares approach to estimate the complex 
amplitude of these P scatterers at each of these L possi- 
ble frequencies. It is also assumed that various attenua- 45 
tions, time delays and angles between the P scatterers 
and I receivers are available via some other means. 
Since, of course, the Dopplers are not known a-priori, 
the method would essentially segment the entire Dop- 
pler range into intervals of delta-f and measure the com- 50 
plex signal amplitude of each of these frequencies and at 
each of the I receivers by some prior known method. 
Essentially, the method obtains the knowledge of the 
spectral contents of the scatterers from the known spec- 


ments from the ship to a plurality of known landmarks. 
The range measurement is done by conventional pulse 
radar. It does not describe any scheme for Doppler 
measurement, which is part of the area of interest ad- 
dressed by this invention. 

U.S. Pat. No. 4,179,696 of Quesinberry contains a 
Kalman estimator which essentially transforms the line- 
of-sight (LOS) measurement data (i.e. range, angle, etc.) 
to obtain the target position as referenced to a stable 
coordinate system. This again is a quite different prob- 
lem. It involves no scheme for a Doppler measurement 
and is essentially concerned with the transformation of 
one set of data (LOS coordinates) into a different set 
(fixed or stable coordinates) as is the case with the Vir- 
not patent mentioned above. 

The patent of Wu (U.S. Pat. No. 4,471,357) describes 
a method of two-dimensional processing of SAR images 
which is an entirely different problem than that of pri- 
mary concern herein. 

The following general commentary may be helpful in 
placing a proper perspective on the present invention as 
will be described hereinafter in relation to least squares 
algorithms and the Kalman filter. The basic idea of 
minimizing the sum of squares of the residuals goes back 
to Gauss in the 19th century who applied it to the study 
of problems in celestial mechanics. Since then, as is 
evidenced by a very large body of published scientific 
research, this generic idea has permeated through many 
scientific, engineering, and other fields. In this regard, 
the celebrated Kalman filter is also a specific implemen- 
tation of the least squares algorithm; but, limited to a 
very structured problem of linear systems with Gauss- 
ian noise. For the more general nonlinear case (as is true 
with the environment of principal concern herein) there 
is no single scheme of a general nature. One possible 
solution, although not always satisfactory, is to use a 
linear approximation and apply the Kalman algorithm 
to this linearized model. A sinusoidal function can 
hardly be considered to be linear except over a very 
small segment. Thus, the above-referenced algorithms, 
although using least squares and Kalman filters, do not 
solve this problem. 

Indeed, the fact that there are so many different tech- 
niques to solve the problem indicates the importance of 
the problem. This, however, also implies that there is no 
single technique superior to all others in all possible 
situations and/or with respect to different criteria such 
as computational complexity, statistical efficiency, etc. 

STATEMENT OF THE INVENTION 


tral contents of the received signals. It does not address 
the problem of how to obtain the spectral contents of 
the received signal. This latter problem is the subject of 
the present invention. Moreover, the framework of the 
scheme is also very computationally intensive as it in- 
volves an exhaustive search in the Doppler domain. It 
also does not take into consideration the time-varying 
Doppler — a clear possibility in a near field situation. 
Thus, in this framework, a discrete scatterer with time- 
varying Doppler would be identified by L complex 
amplitudes that are functions of time rather than identi- 
fying it as a single sinusoidal signal with time-varying 
amplitude, phase and frequency as would be more desir- 
able. 


55 Accordingly, it is the object of this invention to pro- 
vide a novel and elegant least squares scheme for the 
nonlinear estimation problem of estimating the phase, 
frequency, and amplitude of a sinusoid, very rapidly and 
precisely without making any linearizing approxima- 
60 tion. 

Other objects and benefits of the invention will be- 
come apparent from the detailed description which 
follows hereinafter when taken in conjunction with the 
drawing figures which accompany it. 

65 BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 is a simplified functional block diagram of a 
system employing the method of the present invention 
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to estimate the frequency, phase, and amplitude of an 
incoming signal of interest. 

FIG. 2 is a conceptual aid to the understanding of the 
overall algorithm of the present invention. 

DETAILED DESCRIPTION OF THE 
INVENTION 

The methods and apparatus to be described hereinaf- 
ter are based on a new least squares algorithm which 
can be employed to advantage for fast frequency and 
phase acquisition of sinusoids in the presence of noise. 
The algorithm is a special case of more general adaptive 
parameter estimation techniques. The advantage of this 
particular algorithm is its conceptual simplicity, flexibil- 
ity, and application to general situations. For example, 
the frequency to be acquired can be time varying, and 
the noise can be non-Gaussian, non-stationary, and col- 
ored. As the algorithm can be made recursive in the 
number of observations, it is not necessary to have a-pri- 
ori knowledge of the received signal-to-noise ratio or to 
specify the measurement time. This, of course, would be 
required for batch processing techniques such as the 
Fast Fourier Transform (FFT). The algorithm also 
improves the frequency estimate on a resursive basis as 
more and more observations are obtained. When the 
algorithm is applied in real time, it has the extra advan- 
tage that the observations need not be stored. The algo- 
rithm also yields a real time confidence measure as to 
the accuracy of the estimator. 

The inventor herein has authored a paper entitled 
“Fast frequency acquisition via adaptive least-squares 
algorithm” (hereinafter referred to as “the Paper”) 
which describes his early work with respect to his algo- 
rithm. The Paper was published in the Proceedings of 
the International Telemetering Conference, Las Vegas, 
pp. 91-101, Oct., 1986 and later re-published in ex- 
tended and modified form in the IEE PROCEED- 
INGS, Vol. 136, Pt.F, No. 4, Aug. 1989. Therefore, the 
contents of the Paper in its 1986 version only are prior 
art to the present invention which comprises a novel 
improvement to that prior work thereby solving certain 
problems associated therewith. For simplicity, a copy 
of the extended and modified form of the Paper is at- 
tached hereto as Appendix A and the contents thereof 
are incorporated herein by reference. Only those as- 
pects which are an improvement over the disclosure of 
the Paper will be addressed in detail hereinafter. 

The algorithm as disclosed herein differs from the 
present inventor’s earlier published work in the Paper in 
several non-obvious and important ways which will 
now be delineated. 

To keep the truncation errors small, the dimension n 
of the state vector x* in equations (4) and (5) of the 
Paper must satisfy the following bound: 

n = flT 

where 0 is the initial frequency uncertainty in rad/s and 
T is the observation interval. Thus, for example, if 
n=27TXlOO rad/s and T=0.1s, then the dimension n 
must be at least 60. If the observation interval is in- 
creased to Is, n must be greater than 600. Moreover, as 
the computational requirements (in one implementa- 
tion) for processing each signal sample increase as n 2 , 
one would require in the order of 3,600 and 360,000 
operations respectively for these two examples. Obvi- 
ously, this is a very high computational requirement and 
in actual practice one has to simplify the algorithm 
perhaps by some simple ad hoc procedure to keep the 
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computational requirements manageable. However, 
such a possible ad hoc simplification would result in the 
loss of optimality in terms of the estimation errors. 

Related to the foregoing is the fact that the variance 
of the frequency estimation error given by equation (16) 
of the Paper contains a constant K. For low values of 
the state dimension n, the constant K is equal to 4 and 
then the results correspond to the bound for optimum 
variance derived in the literature for the case of high 
signal-to-noise ratio (SNR). For high value of n, how- 
ever, K may be higher. The reason for this is that equa- 
tion (16) is derived when o)d is estimated via equation 
(12) and is based on just two major elements of the 
parameter matrix 6 . For a low SNR case and/or high 
value of n, to keep the optimality, the final estimate of 
o>d should be based on the complete matrix. Thus, a 
procedure is required for optimally combining the ele- 
ments of 0 so as to yield an optimum estimate for c od for 
all conditions. The present invention accomplishes this. 

Since the algorithm is based on minimizing the error 
between the estimated frequency and phase at any par- 
ticular instant in time and the present value of a stan- 
dard against which the incoming signal of interest is 
compared, the computational requirements of the algo- 
rithm as presented in the Paper tended to increase con- 
tinually as the number of terms to be computed also 
continued to increase. This, of course, is not a desirable 
attribute under actual working conditions. That is, the 
present inventor’s prior work provided to be suitable 
for its intended purpose under short term laboratory 
conditions; but, would impose an inordinate computa- 
tional load on actual field equipment employing the 
algorithms thereof. The present invention eliminates 
that problem. 

Finally, it was not obvious or apparent from the dis- 
closure of the Paper how the basic algorithm could be 
extended to the important case of M sinusoids. The 
present invention discloses this aspect for the first time. 

Thus, although the present inventor’s prior publica- 
tions introduced a novel concept, they did not contain a 
simple realizable algorithm except in some limited cases 
as for the case of small values of HT product. The ex- 
tension to the more general problem of estimating sev- 
eral sinusoid was not attempted as the way for such a 
generalization was not clear at the time. The present 
invention provides elegant solutions to all of the above- 
described problems as will now be discussed in detail. 
FIG. 1 shows the manner in which the present inven- 
tion can be implemented. The system 10 inputs the 
signal of interest at 12. The signal 12 is mixed with the 
output of a numerically controlled oscillator (NCO) 14 
and a signal 90° out of phase with the output of the 
NCO 14 at the mixers 16 , 16 ', respectively. As is well 
known, when two frequencies are mixed together, there 
are two outputs from the mixing process, their sum and 
their difference. The lower possible frequency (i.e. the 
difference), of course, also represents the error between 
the two frequencies. Thus, in the system 10 of FIG. 1 , 
the sum frequency signals are removed at 18 and the 
difference frequency signals are sampled at 20 to pro- 
vide quadrature sample signals at 22 and 22 '. The two 
quadrature sample signals 22,22' are then employed at 
24 , 26 , and 28 to compute estimates of the amplitude, 
frequency, and phase of the error signal generally ac- 
cording to the techniques of the Paper. At 30 the NCO 
14 is adjusted in a manner which will tend to drive the 
error towards zero. The error amount is what acts as a 
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real time confidence measure as to the accuracy of the 
estimator as mentioned earlier herein; that is, as the 
computation progresses over multiple samples and the 
estimation closes in on the actual amplitude, frequency 
and phase, the error approaches zero. Thus, the closer 
the error is to zero, the higher the probability that the 
estimate is accurate. A matrix implemented as part of 
the estimation algorithm actually provides an estimate 
of the variance of the estimation error. 

In the algorithm according to the present invention, 
the dimension of the state vector n is determined from 
equation (18) of the Paper and is given by: 


10 


present invention to multiple sinusoids is attached 
hereto as Appendix A and the contents thereof are 
incorporated herein by reference. 

APPENDIX A 

Fast Frequency Acquisition Via Adaptive 
Least-Squares Algorithm 

Consider the problem of estimating an unknown fre- 
quency Wrffrom the measurement y*, z* below: 


n = ^216 — -jj— ft ) 


15 


and thus n is made independent of the observation inter- 
val T. Thus, for 12 = 200 7rrad/s and a typical value of 
P/No24 dB-Hz(251), the required value of n is equal to 20 
8 as against n>60 required before for the case of 
T=0.1s. This reduces the computational requirements 
by a factor of more than 50. For the case of T= Is, the 
improvement factor is 2,500. The important thing to 
realize is that n is no longer a function of T. This is 25 
achieved by segmenting the total observation period 
into To=nf2 _1 s intervals, applying the algorithm to 
successive intervals as described in the Paper, and at the 
end of each interval applying a frequency and phase 
correction to the NCO 14 of FIG. 1 . The overall effect 30 
of this method is that the algorithm effectively seems to 
be operating on a smaller than co^o initial frequency o)di 
between (0, 2To)> an even smaller frequency o)di be- 
tween (0, 3To) and so on without losing any optimality 
in terms of the estimation accuracy. This is because, 35 
according to the present invention, the frequency and 
phase of the NCO 14 are constantly being adjusted 
towards the actual incoming signal of interest such that 
the error is driven towards zero and, correspondingly, 
the number of terms to be computed decreases as the 40 
error also decreases rather than increasing as in the 
prior art approach described in the Paper. This is un- 
doubtedly the major contribution of the present inven- 
tion as it turns what would otherwise have been a very 
limited use algorithm into a most viable approach to 45 
solving the problem at hand on a commercial and prac- 
tical basis. 

The present invention can employ either of two 


yk = A sin (wjtk 
zk ~ A cos (wjtk 


+ 4 > ) + n ik 

>A: = 1, 2, . . . 
+ $) + n qk J 


0) 


Here, {y*, zk } represent the samples of the inphase and 
quadrature components of a. received signal s(t) ob- 
tained by demodulating s(t) by a carrier reference signal 
r(t) and its 90° phase-shifted version respectively, i.e. 


s<f)=/* sin (h>o/ + 4> 0 ) +/;(/) 

rit)—2 sin (K> c r+cJ) c );4> = 4> 0 — <J>c, wj=w 0 — w c 


( 2 ) 


with n,k and n q k denoting the inphase and quadrature 
components of white noise n(t) with variance c r 2 . The 
algorithm can be easily extended to the case where n(t) 
is a colored noise. 

The measurement equations can be written in alterna- 
tive forms as follows: 


yk—A sin (wjtk) cos 4>-M cos (wjtf.) sin <p + »fk 
zfc—A cos (yvdtk) cos 4 >— A sin (wjtk) sin <f> + n g k 


(3) 


or, with a power series expansion for the sine and cosine 
functions, as 




A sin <J> 

A sin 4> ^(cos 4 >)h’j — ^ wf — 


(4) 


A cos <j> — ,4 (sin 4>)wtf — 


2 ! 

A cos <f> 


WJ 2 + 


methods which combine the elements of the matrix 6 so 
as to produce the optimal estimate of od rather than 50 

A cos 4> 
3! 

■ * • 

A sin <£> 

• (n - 1)! 

■*r 

deriving it on the basis of just two elements of 0 via 
equation (12) as was previously done at the cost of 

A sin <J> 
3! 

* Htf 3 . . 

A cos <J> 

■ (« - D! 

■<- 


losing optimality, especially for the cases of large n and 
low SNR. According to one method, by simply adding 
the appropriate terms of 6 one can estimate Sin a>d and 55 
Cos a) d' From a well known theorem in estimation the- 
ory, optimality of 6 also implies the optimality of the 
estimates of Sin coj and Cos o>rf. From these estimates 
a ^ is then uniquely estimated. The alternate method is 
to apply Prony’s method [see N. Mohanty, Random 60 
Signal Estimation and Identification , Van Nostrand 
Reinhold Company, New York, 1986] to optimally 
estimate ft >d from the matrix 0 as shown in the Paper. 
FIG. 2 is included herewith as a conceptual aid in the 
understanding of the overall algorithm of the present 65 
invention. 

A detailed description by the inventor herein relative 
to extending the adaptive least-squares algorithm of the 


t n— 1 
*k 


[:] 


In the above approximation the terms of the order 
(w</t&)Vn! and smaller order have been ignored (assum- 
ing here that wjtjtCn). With obvious definitions, the 
measurement equation can be written in a form ‘linear in 
parameters’: 


Zk-0'xk+nk 


(5) 


In the above, ' denotes transpose, Z'* = [y* zjt], n'^“[n/A: 
x'k denotes the observable state vector [1 1* t* 2 . . . 
t*"” 1 ] and 0 ' is the unknown parameter matrix. A stan- 
dard least-squares algorithm can be applied to estimate 
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the unknown parameter matrix O' from the sequence of 
noisy observations Z*, k= 1, 2, . . . , N. 

The parameter matrix O' can be estimated by either a 
recursive or nonrecursive form. We consider the non- 
recursive form. The estimate of 6 on the basis of mea- 
surement Zk, k= 1,2, . , . , N, denoted On, is given by 


0 0-.'^) 


#A = ^.2^ Xjx'jX^-J ^ 2 ^ xj(x'j6 + n'j)\ x -J ^ 


A simple manipulation of Equation 8 yields the esti- 


mation error 


At A 

0/V ~0 — On as 


where 0<X< 1 is the exponential data weighting factor. 
From §n, the estimates of A, w^and <f> can be obtained. 

The algorithm of Equation 6 requires an inverse of a 15 
symmetric (n X n) matrix once, requiring order n 2 com- 
putations. It may appear that the computation of each 
Xjx'j term requires n 2 computations. However, detailed 
examination shows that only 2n computations are re- 20 
quired. Thus, the total number of computations is equal 
to 6nN + 0(n 2 ). In practice, the matrix inverse can be 
precomputed, thus reducing the data-dependent com- 
putations to only 2nN+(n 2 /2). 

The matrix (lj = \ N Xjx'j\ N -J) in Equation 6 has a 25 

very special structure as can easily be seen by explicit 
computation of the term Xjx'j of the summand. Thus, 


- -If, ) 


As the state vector x; is deterministic and n» is a zero 


VCLIUI Ay 

mean process, $/yhas its^mean equal to zero. The error 
covariance matrix of On can also be evaluated in a 
straightforward manner. Postmultiplying Equation 9 by 
the transpose of @n and taking expected values of both 
sides gives 


= Qlj xjxjkN-j'j 


V t 2 t 2 

v v 0 


j Jj XjE\n'jn,]x'iK^-j> )x ( xpc’^-J ) . 

considering the case of \= 1 and recalling that {n/} is a 
35 white noise sequence 


2n — 2 
’ * 1 V 


Each of the matrices xyx'yand P^ 1 is a Hankel matrix. 
That is, all the elements of each crossdiagonal are the 
same. The structure of Hankel matrix is very similar to 
that of a Toeplitz matrix wherein the elements along the 
various subdiagonals are equal. The well-known fast 
algorithm disclosed in Kumar, “A Fast Algorithm for 
Solving a Toeplitz System of Equations,” IEEE Trans- 
actions, 1985, ASSP-33 (1), pp. 254-267 for the solution 
of a Toeplitz system of equations can be slightly modi- 
fied so as to become applicable to the present problem. 
Thus, Equation 6 can be solved in order n (log 2 n) 2 
computations, resulting in a considerable reduction in 
the requirement for large values of n. 

If the matrix inverse is precomputed then, with the 
well-known fast algorithm referenced above, the solu- 
tion for §n can be obtained in approximately 6n log 2 n 
operations. In the implementations above, it is sufficient 
to store only the first row and column of P or P-l. 

In this case only the measurements {y*} are available 
and the parameter matrix O’ is of dimension nxl. In 
such an implementation, however, there may result an 
ambiguity of tt radians in the phase estimate if the sign 
of wj is also unknown. 

Assuming that the model Equation 5 is exact (the 
dimension n of the parameter matrix is Equation 4 is 
sufficiently high), then the substitution of Equation 5 in 
Equation 6 yields 


/V /V 

£[0A r O' n\ 




E[ |j rij\\ 2 ] = E[njj\ + E[n)j\ + E[n\] = <r 2 (11) 

A simple approximate expression can also be obtained 
for the frequency estimation error when the amplitude 
A is known and uniform sampling is used. The fre- 
quency estimate w^ can be obtained as 

*d,N= ((0A 21 ) 2 + (0^ 2 ) 2 }iA ~ 1 (12) 


When the amplitude A is also unknown, it can be re- 
placed by its estimate given by 

4v={(V) 2 +(V) 2 }i (13) 

In the above expressions, ,0^ denotes the (i,j)th element 
of the parameter matrix 0n< The error variance of these 
elements of interest is given by 


£R<&‘) 2 ] + £I<%) 2 ] t? ) 


65 where K approaches a constant with the increase in the 
numbers of observations N. 

For relatively small errors, the frequency estimation 
error w^at=w £ /— w^has variance approximately 
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^G-vO 


For the case of uniform sampling ty=jT 5 , where T 5 is the 
sampling period. Substituting for t/and letting T=NT 5 
denote the observation period, 


^Ihvv 2 ] = * 


k N(N + l)(2A r 4- 1) 


1 

T? 


1 

A 2 


10 

to their variances, as would be the case for small values 
of W’d- 

To keep the dimension n of the parameter matrix 
small, the following estimation method is proposed. 
Dividing both sides of Equation 16 by fl 2 , where fl 
denotes the maximum possible value of wj, and substi- 
tuting T=n/H, one obtains 


os) 10 = 

a 2 p n 3 


(17) 


Selecting a value of 1/36 for the left hand side of the 
In terms of the unsampled system, if the additive above equation allows one to express the maximum 
noise process has one-sided noise spectral density No, 15 f uncertainty that can be resolved by the algo- 

then ct2=2N 0 /T, Thus, rithm has a function of n. Thus, 


£1^1 = 


6 

P 


~ -T = NT S 


( 16 ) 


n = 


20 


216 


P 

No 


( 18 ) 


where P = A 2 /2 is the received signal power and K has 
value approximately equal to 4 for low values of n. 

We note here that in the derivation of Equation 16, 
the approximation error in Equation 5 has been ignored. 

It is difficult to estimate the error due to such finite 25 
approximation. However, from a few computer simula- 
tions, it appears that for n > w</T = (w/Tj)N, such error 
is small. 

The estimator Equation 12 for w^ from the parameter 
matrix 0, though simple from the computational view- 30 
point, may not provide the minimum possible estimation 
error variance* for w as it ignores the other components 
of ¥. In theory, if the number of terms in the parameter 
matrix is sufficiently large, one may first estimate, by 
direct summation, the terms A cos <f> sin w</, A cos <f> cos 35 
Wtf, A sin <j)sin w* and A sin cf> cos w</ without losing 
any optimality in the estimation, as these terms are re- 
lated linearly to 0. From these estimates one obtains the 
estimates for <f> and wj. Provided that | w<*| <7r/ 2, such 
an estimate for w^is unique. The bound on w^is ensured 4 0 
via a simple normalization described subsequently in 
Section 6 of the paper. It may be note^i that the errors 
associated with various elements of 0 are necessarily 
correlated. 

Alternatively, one may resort to one of the several 45 
possible ad hoc procedures to exploit the additional 
information contained in the estimates w/ 1 , n> 1. In one 
such procedure, in estimating the magnitude of w</, one 
may simply compute the average of the nth roots of 
| w/ 1 1 for n= 1, 2 . . . Perhaps a very systematic proce- 5Q 
dure to obtain w d from % is the application of the least 
square method A to this problem. According to this 
method, with Y,- denoting the estimate of A|wtf|'for 
i=l, 2, . . . , n— 1, the estimate of |w</| is 
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The rationale for selecting the value of 1/36 for 
ElWf/ 2 ]/!! 2 , is as follows. Since the additive noise has a 
Gaussian distribution, one may assume that the fre- 
quency estimation error has a Gaussian distribution 
with its standard deviation denoted by cr(w^). The 
above selection thus ensures that 3o-(w^)<H/2. 

Now, for any given value of n (related to the algorith- 
mic complexity permitted) and the P/No ratio, the value 
of maximum possible initial frequency uncertainty H is 
calculated from the expression 1 8 above. The total esti- 
mation period is then divided into intervals of 
To^nfl^ 1 seconds and the numerically controlled os- 
cillator (NCO) frequency and phase are adjusted at the 
end points of these intervals in the following manner. 

Denoting by w^othe difference frequency (frequency 
of the signal at the phase detector output) at t=0-h» its 
estimate w^o is obtained on the basis of observations y*, 
z&; 1 <k<Ko where Ko=To/Tj is the number of sam- 
ples (assumed to be an integer) in any To seconds inter- 
val and t* — (k— l)Tj. At the end of the first To-seconds 
interval, the NCO frequency is adjusted by w^ such 
that the resulting difference frequency wji in the time 
interval (T, 2T) is simply equal to $do=Wdo=Wd 0 . We 
also note that from the above selection of T 0 , max 
(wrfi)<0.5 max (w^o)- The cumulative phase <f>(t) of the 
signal at the phase detector output is then 


(19) 


a w — 1 A-) 

m = .fj it-i 


n«- 1 AA 

2 TO- 1 

l=\ 


To obtain (including sign) we repeat the above 
procedure with Y, replaced by the estimate of one of A 60 
(sin 4>) w/ A (cos <f>) w^ or an average of these two 
terms. As the estimate of |w^j is independent of the 
phase <f>, it may be more robust and thus it may be 
worthwhile to compute both | w</| and w<*at the cost of 
some redundancy. 65 

It may however be noted that the improvements 
provided by such modifications may not be significant, 
if the expectations of w/>, are relatively small in relation 


<K0 = Hyo7b + wtfi (r - 7o) + 4>o; T < t < IT 

- »>d\t + (h 'do - w dl) 70+4*0 

— »wir + b + 4>o 


Thus, if in addition to the frequency correction a phase 
correction A(f>=w < ^)Tois introduced to the NCO (over 
an interval of To— T 5 to To seconds), then the resulting 
phase is simply equal to w^it+^ofor To<t<2To. Thus, 
the estimation of w<n can be obtained by processing 
observations y*, z*, Ko<k<2Ko using the least-squares 
algorithm with its initial estimate set equal to zero. It 
should be noted that the knowledge about the error 
variance of the previous estimates must also be included 
as a priori information in the least-squares algorithm. 

With 0i denoting the parameter matrix obtained by 
replacing w d by w d\ in the matrix 0 defined earlier in 
Equations 4 and 5, its estimate $1,2*0 obtained on the 
basis of observations y*, Z*; l<k<2Ko is given by 



11 


5 , 165,051 


A 

01 . 2*0 = 


/ 2K 0 ^ 

2 xpc'jhP-^-j 


( 2 * 0 
^=^0+1 


XjZjXlKO-J + 


( 20 ) 


5 


[ji\ jc / t >“ D ’ y 



In Equation 20 Oi^ATois obtained by replacing w^by 0 
(the a priori estimate of w^i) in Wq^ko which denotes the 
estimate of 0\ obtained on the basis of observations y 
z*; l<k<Ko. We note that the maximum product of 
frequency uncertainty and observation period given by I 5 
max w^i. 2To is smaller than ftTo<n and thus the re- 
quired dimension n of the parameter matrix need not be 
increased for the increased time interval 2To- 

The above procedure can be generalized to subse- 
quent To-seconds intervals in a straightforward manner, 20 
while the bound on frequency uncertainty and the 
length of observation period product remains satisfied 
for all t. This is so because, in view of Equation 16 
above, the product E[w 2 ]T 2 in fact decreases with the 
increase in total observation period T. 25 

From the development above it is apparent that the 
estimation error variance is still given by Equation 16 
providing that the peak frequency uncertainty is within 
three times the RMS frequency estimation error. If this 
condition is violated, increased truncation errors may 30 
result. Also, by selecting the data weighting coefficient 
X to be strictly less than 1, the algorithm can track a 
slowly varying input signal frequency wo. In this case, 
for example, wji = \^o+ Aw, where Aw represents the 
change in the input signal frequency wo over the first 35 
To-seconds intervals. Such a variation is taken into ac- 
count in the error covariance matrix in Equations 11 
and 20 by adding an appropriate ‘dynamic noise covari- 
ance matrix’ Q to it at each To-seconds interval. 

In an equivalent implementation, v and hence w^may 40 
be made recursive in each measurement y*, z using a 
standard recursive least-squares algorithm and thereby 
also eliminating the need for a matrix inversion. De- 
tailed calculations of the number of computations per 
sample (a function of n, T, Ko etc.) would dictate the 45 
more specific implementation so as to minimize the 
computational requirements. 

In an alternative approach to keep the value of n fixed 
and small with an increase in the total observation per- 
iod, instead of resetting the frequency reference (mak- 50 
ing a correction in the NCO frequency) the time refer- 
ence is reset to zero. Subsequent observations in Equa- 
tion 1 are now with respect to a different phase-refer- 
ence, say ^ denoted At this stage the observations in 
the second To interval are processed to have the phase 55 
reference $ and are then combined with the first set of 
observations. Equivalently, it is required to postmulti- 
ply the second sum on the right hand side of Equation 
6, obtained for the second subinterval of To, by the 
following matrix: 60 


cos (A<J>) —sin (A<f>) 
sin (A<f>) cos (A<f>) 


where: A<f> = <f> — $and add the result to the correspond- 
ing sum for the first To interval. The first sum on the 
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right hand side of Equation 6 is simply multiplied by a 
factor of 2. This procedure is extended in an appropriate 
manner to subsequent intervals, so as to obtain a final 
estimate for wj and <j> based on the complete set of 
observations. 

Wherefore, having thus described the invention, 
what is claimed is: 

1. A method for accurately and quickly estimating 
signal parameters comprising a signal frequency of a 
signal of interest, said method comprising: 

a) inputting the signal of interest; 

b) generating a reference signal with adjustable refer- 
ence signal parameters comprising a reference fre- 
quency; 

c) mixing the signal of interest with said reference 
signal and a signal 90° out of phase with said refer- 
ence signal to produce a mixed signal and sampling 
said mixed signal to provide an error signal com- 
prising a pair of quadrature sampled signals having 
error signal parameters comprising a difference 
frequency equal to a difference between said signal 
and reference frequencies; 

d) computing estimates from said error signal of said 
error signal parameters employing a least squares 
estimation; and, 

e) adjusting said adjustable reference signal parame- 
ters of said reference signal as a function of esti- 
mates of corresponding ones of said error signal 
parameters in a manner which reduces said error 
signal parameters, whereby said adjustable signal 
parameters approximate said signal parameters of 
said signal of interest. 

2. The method of claim 1 wherein: sampling said 
mixed signal at a sampling rate \/T s during successive 
uniform computation periods To; 

said computing step comprises computing successive 
estimates of said error signal parameters during 
said successive computation periods; and 

said adjusting step comprises adjusting said adjustable 
parameters at the end of each successive computa- 
tion period, whereby each error signal parameter 
generally decreases in successive computation peri- 
ods. 

3. The method of claim 2 wherein: 

said computing step employs a matrix which trans- 
forms between said pair of quadrature signals and a 
state vector which is a function of the time of sam- 
pling by said sampling means, said matrix compris- 
ing coefficients which are functions of said error 
signal parameters obtained by a power series ex- 
pansion, said power series expansion being of order 
n, wherein n is the dimension of said state vector. 

4. The method of claim 3 wherein: 

said computing step comprises producing said matrix 
by least squares estimation and computing said 
estimates from said matrix. 

5. The method of claim 3 wherein: 

said signal parameters comprise the frequency and 
phase of said signal of interest, and wherein n is 
selected to be a function of a maximum uncertainty 
ft in the frequency of said signal of interest and the 
computation period To is selected to be n /ft. 

6. The method of claim 5 wherein said function of ft 
is (ftNoC/P)£; wherein: 

P is the power of said signal of interest, 

7. The method of claim 6 wherein C is approximately 
in the order of 216. 
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8. The method of claim 2 wherein: 
said signal parameters comprise the frequency and 

phase of said signal of interest, and wherein said 
adjusting step comprises adjusting the reference 
frequency of said reference signal by an amount 5 
corresponding to the estimate produced by said 
computing step of the frequency of said error sig- 
nal and adjusting the phase of said reference signal 
by an amount corresponding to the product of said 
estimate of said frequency of said error signal mul- 10 
tiplied by said computation period. 

9. The method of claim 1 wherein: 

said signal of interest comprises plural sinusoids and 
wherein the signal parameters of each sinusoid are 
estimated. 15 

10. A system for accurately and quickly estimating 
signal parameters comprising a signal frequency of a 
signal of interest, said system comprising: 

a) input means for inputting the signal of interest; 

b) a numerically controlled oscillator for generating a 20 
reference signal with adjustable signal parameters 
comprising a reference frequency at an output 
thereof; 

c) mixing means connected to said input means and 
said output of said numerically controlled oscilla- 25 
tor for mixing the signal of interest with said refer- 
ence signal and a signal 90° out of phase with said 
reference signal and sampling means for sampling 
the output of said mixing means to provide an error 
signal comprising a pair of quadrature sampled 30 
signals having error signal parameters comprising a 
difference frequency equal to a difference between 
said signal and reference frequencies; 

d) computation means for computing from said error 
signal estimates of said error signal parameters 35 
employing a least squares estimation; and, 

e) correction means for adjusting said adjustable ref- 
erence signal parameters as a function of estimates 
of corresponding ones of said error signal parame- 
ters in a manner which reduces said error signal 40 
parameters, whereby said adjustable signal parame- 
ters approximate said signal parameters of said 
signal of interest. 

11. The system of claim 10 wherein: 

said sampling means periodically samples the output 45 
of said mixing means at a sampling rate of \/T s 
during successive uniform computation periods T 0 ; 
said computation means computes successive esti- 
mates of said error signal parameters during said 
successive computation periods; and 50 

said correction means adjusts said adjustable parame- 
ters at the end of each successive computation 
period, whereby each error signal parameter gen- 
erally decreases in successive computation periods. 

12. The system of claim 11 wherein: 55 

said computation means employs a matrix which 

transforms between said pair of quadrature signals 
and a state vector which is a function of the time of 
sampling by said sampling means, said matrix com- 
prising coefficients which are functions of said 60 
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error signal parameters obtained by a power series 
expansion, said power series expansion being of 
order n, wherein n is the dimension of said state 
vector. 

13. The system of claim 12 wherein: 

said computation means comprises means for produc- 
ing said matrix by least squares estimation and 
means for computing said estimates from said ma- 
trix. 

14. The system of claim 12 wherein: 

n is selected to be a function of a maximum uncer- 
tainty O in the frequency of said signal of interest 
and the computation period To of said sampling 
means is selected to be n/H. 

15. The system of claim 14 wherein said function of fl 
is (ftNoC/P)i, wherein: 

P is the power of said signal of interest, 

No is the noise spectral density, and 

C is a selected constant numerical value. 

16. The system of claim 15 wherein C is approxi- 
mately in the order of 216. 

17. The system of claim 11 wherein: 

said signal parameters comprise the frequency and 
phase of said signal of interest, and wherein said 
correction means adjusts the frequency of said 
reference signal by an amount corresponding to the 
estimate produced by said computation means of 
the frequency of said error signal and adjusts the 
phase of said reference signal by an amount corre- 
sponding to the product of said estimate of said 
frequency of said error signal multiplied by the 
sampling period of said sampling means. 

18. The system of claim 10 wherein: 

said signal of interest comprises plural sinusoids and 
wherein the signal parameters of each sinusoid are 
estimated. 

19. The system of claim 10 wherein: 

said computation means provides said error signal 
parameters as a real time confidence measure as to 
accuracy of said estimates whereby the closer said 
error signal parameters are to zero, the higher the 
probability that said estimates are accurate. 

20. The system of claim 10 wherein: 

said computation means includes a matrix as part of 
an estimation algorithm implemented thereby 
which provides an estimate of the variance of an 
estimation error. 

21. The system of claim 1 wherein: 

said computing step provides said error signal param- 
eters as a real time confidence measure as to accu- 
racy of said estimates whereby the closer said error 
signal parameters are to zero, the higher the proba- 
bility that said estimates are accurate. 

22. The system of claim 1 wherein: 

said computing step employs a matrix as part of an 
estimation algorithm implemented thereby which 
provides an estimate of the variance of an estima- 
tion error. 

***** 
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